arXiv:1506.08817vl [astro-ph.HE] 29Jun2015 


Limits on anisotropy in the nanohertz stochastic gravitational-wave background 
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The paucity of observed supermassive black hole binaries (SMBHBs) may imply that the gravitational 
wave background (GWB) from this population is anisotropic, rendering existing analyses sub-optimal. 

We present the first constraints on the angular distribution of a nanohertz stochastic GWB from circular, 
inspiral-driven SMBHBs using the 2015 European Pulsar Timing Array data (Tj. Our analysis of the GWB 
in the ~ 2 — 90 nHz band shows consistency with isotropy, with the strain amplitude in 1 > 0 spherical 
harmonic multipoles < 40% of the monopole value. We expect that these more general techniques will 
become standard tools to probe the angular distribution of source populations. 

PACS numbers: 04.80.Nn, 04.25.dg, 95.85.Sz, 97.80.-d 97.60.Gb 04.30.-w 


Introduction.- Pulsar Timing Arrays (PTAs) are cur¬ 
rently being used to search for, and to eventually charac¬ 
terize, the nanohertz stochastic gravitational-wave back¬ 
ground (SGWB) by looking for correlated deviations in the 
pulse times of arrival (TOAs) of multiple radio millisec¬ 
ond pulsars distributed across the sky. The SGWB in the 
nanohertz regime is thought to be generated by the inco¬ 
herent superposition of a large number of weak and unre¬ 
solved GW sources, including supermassive black hole bi¬ 
naries (SMBHBs) lIlHSl, decaying cosmic-string networks 
or primordial GWs ifT^ IT^ . Previous analyses 
have assumed background isotropy, which emerges as a 
special case from the more general anisotropy framework 
presented here. Although GWs have not yet been di¬ 
rectly detected, limits on the angular power distribution 
of a nanohertz SGWB may constrain the distribution of 
low redshift structure lITSl . the location of several partic¬ 


ularly bright nearby sources dominating the signal strain 
budget flhl im . and open a new avenue to explore the 
population characteristics of SMBHBs. Moreover, if a 
significant fraction of SMBHBs stall rather than merge, 
or are rapidly driven to merger via strong couplings to 
the galactic nuclear environment, then we may expect a 
depleted nanohertz GW signal dominated by only a few 
bright sources ifTSlI . As such, the tools implemented here 
may provide new and novel insights into the final-parsec 
problem (see e.g. Ref. 091). This research is a result of the 
common effort to directly detect gravitational waves using 
pulsar timing, known as the European Pulsar Timing Array 
lEPTA.IMl. 

Limits on the SGWB are usually reported in terms of 
the characteristic-strain spectrum hc{f) of a background 
which is composed of purely GW-driven, circular, inspiral- 
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ing SMBHBs which obeys a simple power-law: 

Kif) = Ahif, (1) 

where is the strain amplitude reported at a reference 
frequency / = yr“^ ll2Tl . The correlations induced hy 
a SGWB in pulsar TOAs can he understood hy consid¬ 
ering a perturhation to the space-time metric along the 
Earth-pulsar line of sight causing a change in the perceived 
rotational-frequency of the pulsar ll22l - E5l . The fractional 
frequency shift 5v{t)/vQ of a signal from a pulsar at rest 
frequency t'g is the difference in the metric perturhation 
at the Solar System harycenter (SSB), and at the pulsar. 
This frequency shift is integrated over time to give the in¬ 
duced timing residuals, r{t) = 6iy{t') / v^dt' , which are 
cross-correlated between pulsars in an effort to boost the 
detection probability of GW signals at the Earth. The ex¬ 
pectation value of the cross-correlated timing residuals be¬ 
tween pulsars a and b is proportional to the overlap reduc¬ 
tion function (ORE, Fab)- a dimensionless function which 
quantifies the response of a pair of pulsars to the stochastic 
GW background ll2^ |271 . In this Letter, we use analyti¬ 
cally computed anisotropic OREs Il28ll29l and recently de¬ 
veloped Bayesian techniques ll^ to constrain the angular 
power distribution of the SGWB. 

Fitting a pulsar timing model- The average pulse pro¬ 
files of millisecond pulsars are remarkably stable and re¬ 
producible. This stability permits high-precision timing, 
which is crucial to GW searches: the minimum detectable 
GW strain is he oc 10“^®(iTrms/100 ns)(T/10 yr)“\ 
where CTi-ms is the root-mean-square of the pulsar timing 
residuals, and T is the total observation timespan Q. 
Therefore high timing precision and long-term observa¬ 
tions are required to distinguish the GW signal from noise, 
as well as boost the signal-to-noise ratio (SNR) in a search. 

Pulsar observations lead to a catalogue of TOAs which 
can then be analyzed to search for GWs. A timing 
model describing all deterministic contributions to a pul¬ 
sar’s TOAs (rotational frequency, spindown rate, disper¬ 
sion measure (DM), etc.) is iteratively fit with the analysis 
package, TEMP02 lISTlI^ . The difference between the 
measured TOAs and the refined timing model prediction is 
the post-fit timing residual, which constitutes the input data 
in our GW analysis. 

GW analysis pipeline- We use the signal modeling tech¬ 
niques described in Ref. ll^ (from hereon L15). The pos¬ 
terior probability of the model parameters, 0, given the 
concatenated post-fit timing residuals from all pulsars, 5t, 
is a multivariate Gaussian: 

exp 

-^detl2,(QTCG)]- ® 

where p(0) is the prior probability distribution of model 
parameters, and projecting all quantities with the matrix 
G marginalizes this posterior probability over all timing 
model parameters (see Ref. [(34l ). The covariance of the 


pre-fit timing residuals is defined as C = Cred + N 
where Cred includes the SGWB, intrinsic pulsar red noise, 
and DM variation components, while N denotes all white- 
noise components. This red covariance is decomposed 
in terms of a low-rank approximation such that Cred = 
Fwhere F is a block-diagonal matrix of Eourier ba¬ 
sis vectors, and is a spectral covariance matrix ll35] - [y7l . 
Intrinsic red-noise and SGWB are expanded in the same 
Eourier basis, while the DM-variation signal is expanded 
in basis-functions which differ only by an extra multiplica¬ 
tive factor of oc 1 where Vq is the observing frequency 
of the pulses. The matrix N is diagonal, with entries given 
by the squared TOA errors which have been corrected by 
previous single pulsar analyses ll33ll . We apply a multi¬ 
plicative factor to all error bars of a given pulsar (referred 
to as the GEEAC parameter) which is searched over here. 

The matrix p has band-diagonal structure, since Eourier 
modes between different pulsars may be correlated due to 
the presence of a SGWB or correlated noise. Therefore 

F “t“ Pid (3) 

where i and j index the discretely sampled signal or 
noise frequencies in our analysis of pulsar TOAs; p = 
/ic(/)^/(127r^/^Tinax) IS the power spectrum of the 
SGWB, with Tinax equal to the timing baseline of the PTA; 
e is the spectrum of a completely correlated red-noise pro¬ 
cess which may result from modeling inaccuracies due to 
drifts in the observatory and global time standards; p is the 
spectrum of a common, but uncorrelated, red-noise process 
which may originate from common physical processes in¬ 
side the neutron stars, see e.g. Refs. |[38]|39l; and Ka is the 
individual red-noise and DM-variation spectrum for pul¬ 
sar a. All these spectra are modeled with power laws, 
(A7l27r2rj„ax)(/n/yr"^)^“"^ yr^ where A = A^ for 
the SGWB; a is a spectral index which equals —2/3 for the 
SGWB; and /„ are the n frequencies at which we sample 
the spectra of red-noise processes, where in this analysis 
n = 50. 

The ORE, Fab, is the average of the overlap of the 
pulsars’ antenna response functions, F^{Cl) (see e.g. 
Refs. 1281 |30l), over GW propagation directions O, and 
weighted by the SGWB angular power distribution, F’(O): 

rab = ^ii +Sab)[ dOP(o) (4) 

Svr Js2 ^ 

where q labels the {-f, x} GW polarization. An ex¬ 
cess in F’(O) in a particular region of the sky may in¬ 
dicate a particularly bright single source, or a hotspot 
of several sources l28l [30l [40j |4T1. In the following 
we decompose the SGWB angular distribution such that 
P{n) = with normalisation 

Jg 2 P{Cl)di). = dvr, where are the real spherical har¬ 
monics. Inserting this decomposition into Eq. Q and pro¬ 
ceeding as in Ref. ll28ll . we expand Fq;, into a sum over 
anisotropic OREs, F|^^\ with associated weights, C/^. to 
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be constrained by the analysis and which characterize the 
SGWB angular power distribution. We note that the lead¬ 
ing function in this expansion, CoorgQ^\ corresponds to the 
ORF applicable to the monopole moment of P{Cl) (also 
known as the Hellings and Downs curve ||42l). Hence, cur¬ 
rent analysis strategies for isotropic SGWBs emerge from 
our fully general anisotropy framework as a special case. 

Each pulsar has 5 stochastic parameters to be con¬ 
strained in a Bayesian analysis: intrinsic red noise (^4, a), 
DM-variation {A, a), and a GEFAC parameter. The fully 
correlated red-noise component, e, will contribute 2 power- 
law parameters, as will the common, uncorrelated process, 
rj. The spectrum of the SGWB is modeled with a fixed 
slope of —2/3 and an amplitude, to be constrained, 
while (max > 0 analyses will include [((max + 1)^ ~ 1] 
additional parameters. The priors on parameters are: 
logigA e (7[-20,-10], a G [/[-2.0,1.5], GEFAC G 
U[0.1, 10.0]. The Cim coefficients are constrained by a 
prior requiring the implied distribution of GW power to be 
positive at all sky-locations |[30l, called the physical prior. 
The prior on Ah is treated separately from other red-noise 
components, and is uniform in the range [10“^°, 10“^°]. 
Applying a uniform prior on Ah with logarithmic priors 
on the amplitudes of all other red components will provide 
the most conservative upper-limits on the strain-spectrum 
amplitude of the SGWB. 

Results.- We parametrize the angular distribution of the 
SGWB down to the angular resolution of the PTA. The 
most anisotropic SGWB signal is one dominated by a sin¬ 
gle source. Hence the angular resolution, and thus (max. 
is a function of the number of pulsars, which sig¬ 

nificantly contribute to a single-source detection, and the 
SNR of that detection |i43l . Sesana and Vecchio ll^ 
find that the angular resolution of a PTA for a resolvable 
GW source is AG oc 50 (50/Apsr)^^^ (10/SNR)^ deg^, 
and this resolution sets an upper bound on ( via ( = 
180/0, where 9 = s /AG Il29l . We analyse a subset 
of the six best pulsars in the EPTA U) which encapsu¬ 
late ~ 95% of the full-array SNR in simulated continu¬ 
ous GW searches gl: PSRs J0613-0200, J1012-h5307, 
J1600-3053, J1713-f0747, J1744-1134, J1909-3744, 
where Tmax = 17.7 years and the GW frequencies 
with which we characterize red-noise components are 
G [l/Tmax = 1-79 nHz, 50/Tmax = 89.7 nHz]. Hence, 
in our array subset (max ^ 4. Carrying out searches 
with the noise characteristics of each pulsar fixed, we find 
the upper limits on the strain amplitude remain consistent 
whether we analyze this six-pulsar subset or the full ar¬ 
ray. Including more pulsars of comparably high timing 
quality would contribute a larger number of pulsar-pairs 
[Apairs = Apsi.(Aps,. — l)/2], which would serve to in¬ 
crease the SNR and resolving power ((max) of any search 
for anisotropy. This comes at the cost of longer likeli¬ 
hood evaluation times, making the systematic study pre¬ 
sented here currently intractable. Our goal is to provide 


the first constraints on anisotropy in the SGWB via a sys¬ 
tematic study with current techniques - we do so with the 
15 distinct pulsar-pairings afforded by a six-pulsar array. 
All analysis is performed with parallel-tempering Markov 
chain Monte Carlo (MCMC) analysis. 

The 95% upper limits on Ah from our analyses are 
shown in Table [I] Firstly, we perform searches with a sin¬ 
gle set of anisotropy coefficients Cim across the entire band, 
which we call the {All-band anisotropy parametrization, cf. 
Table [^. We also perform frequency-dependent searches 
by parametrizing each frequency with independent Cim{f) 
coefficients. We split our band into 5 equal sub-bands 
(A/ = 17.9 nHz), and independent Cim{f) coeffi¬ 

cients constrained in each, called the Frequency-dependent 
anisotropy parametrization (i). Finally, motivated by the 
results of F15 - where most of the SGWB constraints 
were found to come from the lowest three frequencies- 
we apply independent C/m(/) coefficients to the lowest 
four frequencies in our analysis (/ = [1, 2, 3,4]/rmax = 
[1.79,3.59,5.38,7.18] nHz), with the remainder of the 
band(/ = [5,...,50]/rmax = [8.98,..., 89.7] nHz) 
parametrized by a single set of coefficients. This is re¬ 
ported as Frequency-dependent anisotropy parametriza¬ 
tion (ii). The recovered upper limit does not deteriorate 
through the increased number of parameters in our higher 
multipole searches. The monopole upper limits do not pre¬ 
cisely match those found in F15 due to several variations 
in the analysis specifics, namely: (i) our prior on the am¬ 
plitude of red-noise components is uniform in logarithm of 
the amplitudes which provides a more conservative upper 
limit on the SGWB strain amplitude; (ii) we do not con¬ 
sider Solar System ephemeris errors in our correlated noise 
modeling; and (in) we employed a pure time-domain like¬ 
lihood in the initial single pulsar analysis to correct the 
TOA errors in each pulsar. Hence, our monopole upper 
limits are higher than in F15 by ~ 1 x 10“^®. However, 
moving beyond the first analysis presented here, our more 
general anisotropy framework can be easily incorporated 
into all existing and planned pipelines to become a standard 
toolset, since it recovers the isotropic SGWB constraints as 
a special case. The upper limits on the strain-amplitude 
in each anisotropic multipole of the search are shown in 
the left panel of Fig. [^ where the constraints are entirely 
dominated by the restrictions imposed on the by the 
physical prior. Our data are not informative enough to up¬ 
date the prior knowledge we have about the anisotropy of 
the GW sky. 

Rather than impose a specific decomposition of the 
SGWB sky during sampling, we can recover the cross¬ 
correlation values between pulsar pairs and map these to a 
chosen basis in post-processing. We perform a Bayesian 
search for the distinct elements of the Cholesky factor 
of the residual cross-correlation matrix, which ensures 
positive-definiteness of the final matrix ll35l |45]| . After 
sampling we define a mapping between the coefficients of 
the ORF in a particular basis, c, and the cross-correlation 
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FIG. 1: 95% upper limits on the strain amplitude, where (7; = J2m=-i lGml^/(2^ + !)■ Left- all-band anisotropy parametrization 
and frequency-dependent parametrization (ii). The right axis is the ratio of the upper limit to the monopole. The inset figure shows 
95% upper limits on which are marginalized over the strain amplitude for the all-band anisotropy parametrization and a 

constant likelihood analysis. Our limits reflect the constraints of the physical prior. Right: all-band anisotropy parametrization, where 
the cim values are obtained by mapping cross-correlation values to the spherical harmonic basis, without physical prior rejection. 


values, r, such that F = He. A single row of the ma¬ 
trix H will have entries corresponding to the ORF be¬ 
tween pulsars a and b evaluated for all basis terms. In 
the spherical-harmonic basis, such a row would consist 

of and for a pixel basis this is 

Having recovered posterior sam- 

\ iil ii2 ^ 

pies of the vector F, we map these to samples of c via 
c = H^F, where H+ corresponds to the Moore-Penrose 
pseudo-inverse of the matrix H ll46l WT\ . The results 
for mappings to the spherical-harmonic basis with vary¬ 
ing ^max are shown in Fig. [fright). The data support 
such strong anisotropy signatures in this model because 
the joint-posterior in the cross-correlation values are con- 
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FIG. 2: 95% upper limits on the GW strain amplitude in each 
pixel. These limits are obtained by mapping from the Bayesian 
MCMC-sampled cross-correlation values to a pixelated ORF ba¬ 
sis (Apix = 12288). White stars show the pulsar locations. 


sistent with essentially the entire range of [—1,1], which 
when mapped to a spherical-harmonic ORF-basis leads to 
large Cim values. There is nothing to penalize these large 
anisotropy coefficients, which lead to highly anisotropic 
(and possibly negative) GW power distributions and would 
otherwise be restricted by the physical prior. This supports 
to our claim that the constraints in Fig. [T] (left) are prior- 
dominated. 

We also map our recovered cross-correlation samples to 
a pixel basis with 12288 equal-area pixels on the sky. We 
supplement our mapping with the additional normalization 
constraint that P{^)dCl ~ c{Q.i)ACli = dvr. 

The resulting SGWB power in each pixel is marginalized 
over all other pixels and truncated to obtain the positive 
ID-marginalised power PDF before it is integrated over to 
obtain the upper limit on the strain-amplitude in that pixel. 
The result is shown in Fig.[^ where we see the distinc¬ 
tive overlapping antenna patterns of the pulsars mapping 
out the sensitivity of the PTA to the background strain- 
amplitude. The constraints on from each pixel are quite 
poor, and in some cases are more than an order of magni¬ 
tude worse than the all-sky upper limit. As we decrease 
the resolution of the pixelation the constraints in each pixel 
become tighter, until we reach the limit of one pixel, which 
recovers the usual all-sky upper-limit. Figure can also 
help to explain the results in the right panel of Fig.[T] where 
we see that the distribution of pulsars in our array leads to 
the sub-optimal overlapping of the antenna response func¬ 
tions, which in turn causes insensitivities around the 4 clus¬ 
tered pulsars and on large angular scales. Hence, we will 
lack sensitivity to large angular scale anisotropy (L ~ 1), 
which is reflected in the right panel of Fig.[2 Moreover, 
this sensitivity map illustrates the importance of timing 
pulsars from all over the sky to ensure a more uniform sen¬ 
sitivity to GW strain, which will be possible through inter- 
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^max 

Ah, all-band c/^ 

C/m — C/?7i(/)^ ^ 

Ah', Clm = CZm(/)^”^ 

0 

3.94 X 10“^^ 

N/A 

N/A 

1 

4.09 X 10“^^ 

4.06 X 10“^’" 

4.06 X 10“^^ 

2 

4.06 X 10“^® 

4.07 X 10“^’" 

4.02 X 10“^® 

3 

4.06 X 10“^® 

3.98 X 10“^’" 

4.01 X 10“^® 

4 

4.03 X 10“^^ 

3.95 X 10"^'’ 

3.99 X 10“^^ 


TABLE I: 95% upper limits on SGWB strain amplitude A^. The 
first column is the all-band anisotropy parametrization, the sec¬ 
ond and third correspond to the frequency-dependent anisotropy 
parameterizations (i) and (ii) respectively, described in the text. 


national collaborations such as the International PTAP48I. 

Conclusions - Our analyses suggest that this dataset is 
not informative enough to update our prior knowledge of 
the angular distribution of the nanohertz SGWB. Using 
a prior which enforces a positive SGWB distribution, we 
find that the 95% upper limit on the strain amplitude in 
multipoles of the background distribution with / > 0 is 
< 40% of the monopole strain. No evolution of these 
upper-limits as a function of GW frequency is found since 
the constraints are a reflection of the prior. Addition¬ 
ally, we can recover the joint posterior distribution of the 
cross-correlation values between pulsar pairings, and sub¬ 
sequently map these to a spherical-harmonic or pixel ORF- 
basis. With the only constraint being positive-definiteness 
of the cross-correlation matrix, the strain-amplitude in I > 
0 multipoles is < 400% of the monopole value. The strain- 
amplitude upper-limits as a function of location on the sky 
reflect the overlapping antenna pattern behavior of the full 
PTA, where the limits can often be more than an order of 
magnitude worse than the all-sky limit. A full description 
of all techniques employed here, and their efficacy, will be 
provided in a follow-up mefhods paper. 

Forfhcoming advanced radio insfrumenfs such as fhe 
Five-hundred-mefre Aperfure Spherical Radio Telescope 
[FAST, |49l l50l . MeerKAT lISTl . and fhe Square Kilome- 
fre Array [SKA,|52l will enhance fhe defection and infer¬ 
ence prospecfs for anisofropic GW skies by defecting large 
numbers of millisecond pulsars and timing them to un¬ 
precedented precision. Upcoming studies will investigate 
how we can combine galaxy catalogues with frequency- 
dependent maps of the nanohertz GW sky to probe whether 
the strain budget is being dominated by a few bright nearby 
sources, or is more diffuse. We hope that the work pre¬ 
sented here, together with these future studies, will provide 
important insights into the demographics, evolution, and 
assembly of SMBHBs not accessible by any other means. 
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